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Summary 


During this reporting period, theoretical work on the secondary atomization 
process was continued and the experimental apparatus was improved. A one-dimensional 
model of a rocket combustor, incorporating multiple droplet size classes, slurry 
combustion, secondary atomization, radiation heat transfer, and two-phase slip between 
slurry droplets and the gas flow has been derived and a computer code has been written to 
implement this model. The STANJAN chemical equilibrium solver has been coupled with 
this code to yield gas temperature, density, and composition as functions of axial location. 
Preliminary results indicate that the model is performing correctly, given current model 
assumptions 

Radiation heat transfer in the combustion chamber is treated as an optically-thick 
participating media problem requiring a solution of the radiative transfer equation. A 
cylindrical P, approximation has been employed to yield an analytical expression for 
chamber- wall heat flux at each axial location. 

The code was exercised to determine the effects of secondary atomization 
intensity, defined as the number of secondary drops produced per initial drop, on chamber 
burnout distance and final A1 ? 0 3 agglomerate diameter. These results indicate that only 
weak secondary atomization is required to significantly reduce these two parameters. 
Stronger atomization intensities were found to yield decreasing marginal benefits. 

The experimental apparatus was improved to reduce building vibration effects on 
the optical system alignment. This was accomplished by mounting the burner and the 
transmitting/receiving optics on a single frame supported by vibration-isolation legs. 
Calibration and shakedown tests indicate that vibration problems have been eliminated and 
that the system is performing correctly. 
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OVERALL PROGRAM OBJECTIVES 


The overall objective of this project is to provide an increased understanding of the 
secondary atomization characteristics of Al/RP-1 slurry propellants. Specific objectives 
are as follows: 

1 . Develop an experimental system to measure the size and velocity of burning 
10-100 pm slurry droplets and to check for the presence of burning aluminum 
in these same droplets. 

2. Use this experimental system to determine the ignition and secondary 
atomization characteristics of various slurry formulations. 

3. From the experimental data develop an understanding of the role surfactants, 
gellants, and ultimate particle play in the secondary atomization process. 

4. Develop analytical models of droplet ignition and secondary atomization. 

5. Apply these models in a 1-D rocket model to determine the effects of 
secondary atomization on engine performance. 
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PROGRESS DURING REPORTING PERIOD 


THEORETICAL EFFORTS 
One Dimensional Rocket Model 

Recent performance analyses of slurry fueled rockets 1 ' 5 predict that A1 slurry 
propellants may increase maximum payload over neat liquid systems. However, these 
studies neglect additional radiation and two-phase flow losses resulting from solid 
combustion products and increases in droplet combustion times due to solid agglomerate 
combustion. Accounting for these three factors is critical in accurately predicting the 
effects of using metallized propellants as rocket fuels. 

Our previous research involving slurry droplet combustion 6 ' 7 focused on the 
mechanism of secondary atomization as a possible means of reducing combustion times 
and two-phase flow losses. However, this research was devoted to understanding the 
droplet combustion and secondary atomization processes themselves and did not involve 
the effects of secondary atomization on engine performance. 

Consequently, a one-dimensional engine model is being derived to provide a 
preliminary evaluation of the effects of secondary atomization, two-phase flow losses, and 
radiation heat transfer on engine performance. Coding a new model provides some 
benefits over using a standard code such as TDK. First, the new model will permit the 
inclusion of secondary atomization effects, and second, the model will allow a clearer, if 
simpler, understanding of the physics involved in the problem. 

Over the past six months, a numerical code incorporating the following effects has 
been written to model a rocket combustion chamber. 

• Multiple droplet size classes 

• Gas-phase chemical equilibrium 
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• Two-phase flow losses 

• Hydrocarbon evaporation and burnout 

• A1 combustion with surface condensation of AI 2 O 3 

• Simple secondary atomization 

Chamber dimensions and propellant flow rates have been taken from Galecki's work 8 to 
provide realistic model inputs. 

Model Description 

The propellants used in the engine model are a JP-10/A1 slurry and a preheated 
gaseous O 2 oxidizer. JP-10, a pure hydrocarbon (C 10 H 16 ), was chosen as the slurry 
hydrocarbon component instead of RP-1 to avoid the complexity of modeling multi- 
component droplet combustion. 

The combustor flow is modeled using a single product-phase containing both gases 
and small A1 2 0 3 fume particles, and three additional flow phases for each droplet size 
class: a liquid hydrocarbon phase, an aluminum phase, and a phase containing large A1 2 0 3 
agglomerates. The AI 2 O 3 must be separated into two flow phases because of the two 
oxidation mechanisms inherent in aluminum combustion . 9 ' 11 In the first oxidation 
mechanism, large A1 2 0 3 agglomerates are formed through droplet surface 
condensation/oxidation, and in the second, very small A1 2 0 3 fume particles are produced 
through vapor-phase oxidation. Since the A1 2 0 3 agglomerates are much larger than the 
fume particles and are attached to the surface of the aluminum droplets it is necessary to 
model the A1 2 0 3 agglomerates separately from the fume particles. Including the fume 
particles in the gas-phase flow requires the assumptions of no temperature or velocity slip 
between the fume particles and the gases, but greatly simplifies the problem solution. 

As an A1 droplet bums, it continually produces A1 2 0 3 fume particles at the oxide 
boiling point which subsequently equilibrate with the gas-phase temperature and velocity. 
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This process results in a range of fume particle temperatures and velocities as new 
particles are created and older particles continue to equilibrate with the gases. If the 
temperature and velocity slip were not neglected, many additional flow phases would be 
required to accurately model this range of fume particle temperatures and velocities 

Verification of No Slip Assumptions 

A1 2 0 3 fume particles formed during aluminum combustion in a rocket environment 
are typically on the order of one micron in size 9-12 Parry and Brewster 12 measured a fume 
particle mean diameter, D 32 , of 0.97 pm. Using the model for particle acceleration 
developed later in this report, it was found that these fume particles rapidly accelerate to 
the gas-phase velocity thus justifying the no-slip velocity condition. 

The no-slip temperature condition requires that the convective heat transfer 
resistance be much less than the radiation heat transfer resistance for a minimal 
temperature difference between the steady-state particle and gas temperatures. The 
combined radiation/convection heat transfer rate must also be large enough to allow 
particle temperature, T p , to rapidly adjust to a changing gas temperature, T gas The 
following is an analysis verifying the no slip temperature condition, which is defined as: 

It -t|<ooit (i) 

At equilibrium, the energy radiated from an A1 2 0 3 particle to the surrounding 
medium equals that convected to the particle from the gas flow. However, the solution of 
a participating medium radiation problem is complex with radiant intensity dependent on 
particle location in the chamber. As a conservative approximation, medium participation 
is neglected and the fiime particles are assumed to radiate directly to the chamber walls 
Using this approximation, the particle energy balance can be expressed as 
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Ah(T„-T p ) = Aeo(i;-i:) 


( 2 ) 


where A is the particle surface area; h, the convective heat transfer coefficient; e, the 
particle emissivity; and o, the Stefan-Boltzmann constant. The wall temperature in the 
radiation term, T w , is an unknown parameter dependent on chamber cooling design but is 
typically much less than the oxide fume temperature and can therefore be neglected 
Substituting for h using the definition of the Nusselt number 


Nu = 



( 3 ) 


where D p is particle diameter, and for T p using the no-slip definition, yields 


-0.99TJ = ea(0.99Tj 4 . (4) 

“P 

Equation 4 can be simplified to yield the maximum gas temperature, above which, the no- 
slip temperature condition is invalid. This yields 


gas, max 


0.0104 


Nuk^ 

D p ec 


( 5 ) 


The fume particles can be modeled as spheres in a stagnant environment by 
employing the no-slip velocity condition, yielding a Nusselt number of 2.0. The fume 
particles can be approximated as uniform 0.97 pm diameter spheres and the gas 


7 



conductivity, k gas , as that of oxygen (k = 0.215 W/m-K at 4000 K). 13 Konopka, Reed, 
and Calia 14 report that impurities can influence the infrared optical properties of A1 2 0 3 
particles produced by solid rocket motors and that these optical properties can vary from 
engine to engine. Because of this uncertainty, a particle emissivity of l 0 is used as an 
extreme value 

Using the above assumptions in Eq. 5, a maximum allowable gas temperature of 
approximately 4330 K is calculated. Currently, the chamber model is adiabatic, predicting 
a worst-case maximum temperature of 3872 K. Accounting for flow heat losses, a 
participating medium and using a fume particle emissivity less than 1 .0 should increase this 
difference even more, clearly indicating that the steady-state fiime particle temperature is 
close to the gas temperature for a constant gas temperature situation 

The fume particle temperature must also respond quickly to changes in the gas 
temperature for the no-slip temperature condition to apply. This second requirement is 
verified as follows: 

Approximating the fume particles as uniform temperature spheres requires that the 
Biot-number criterion be satisfied. 

hL D„ 

Bi = — < 0. 1 where L = — - (6) 

k„_, 6 


Assuming a spherical particle moving at the gas-phase velocity, the Biot number reduces 
to 


Bi = — - kg “ 


3 k, 


(?) 
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where k is thermal conductivity. The minimum thermal conductivity of A1 2 0 3 (-5,66 W/m 
•K at 1500 K) 15 and the maximum gas thermal conductivity (~ 0.22 W/m-K for 0 2 at 
4000 K) 13 were used in estimating the maximum Biot number. This results in Bi « 0.01, 
indicating that A1 2 C>3 fume particles can be modeled using a lumped-temperature analysis. 

Based on the above lumped-temperature approximation, the change in sensible 
energy of a particle can be determined through the unsteady energy equation, 

"’S^ L =-A[l'(T p -T^) + eo(i;-T;)] , (8) 

where m is the mass of an A1 2 0 3 fume particle; c p , the specific heat of A1 2 0 3 ; and A, the 
particle surface area Given an initial particle temperature, T p i nit =0.99T gas jnjt , and a new 
gas temperature, T gas new , the time required for the particle temperature to reach 
0.99T gas ncw can be determined by numerically integrating Eq. 8. If this particle response 
time is less than the time required for combustion to raise the chamber gas temperature 
from T gas ^ to T new , then the no-slip temperature condition is valid for all times. Using 
the following approximate values, p p = 3970 kg/m 3 , D p = 0.97 pm, c p = 1225 J/kg-K, 

Nu = 2.0, s = 1.0, and k gas = 0.2 W/m-K, yields a particle response time of 5.5 ps for a 
sudden jump in gas temperature from 3000 K to 3300 K. In comparison, model 
calculations predict that 250 ps are required for chamber gas temperature to rise from 
3000 K to 3300 K indicating that particle response time is short enough for the no-slip 
temperature condition to apply. 

Governing System Equations 

Mass Conservation : For M droplet size classes, the following equation can be written for 
system mass conservation assuming steady-state conditions and that no mass is added to 
or removed from the chamber except at the injector face and the chamber exit: 
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( 9 ) 


drt V _ yf dm LH | drn^ | "j 

dx dx dx dx J. 

In Eq. 9, x is axial position in the combustor; j, a particular droplet size class; and rh^, 

the gas mass flux. The liquid hydrocarbon, aluminum, and A1 2 0 3 agglomerate mass fluxes 
for a particular droplet size class are given as m LH , m^ , and m Alj0j , respectively. The 

terms on the right hand side of Eq 9 are found from a hydrocarbon droplet gasification 
model and an aluminum droplet combustion model. The development of all three terms is 
similar and is presented below. 

The change in mass of a single slurry droplet in a size class, j, during a time 
interval dt is: 


m 


change , j 


dm 


drop.j 


dt 


dt 


( 10 ) 


where m change j represents the hydrocarbon, aluminum, or A1 2 0 3 mass change. The value of 
dm^opj/dt for the liquid hydrocarbon is found from the hydrocarbon gasification model 
while dn ifr o p j/dt for the aluminum and A1 2 0 3 are both determined from the aluminum 
combustion model Relating dt to dx through the chain rule and the velocity relationship 


u 


p.j 


dx 

dt~ 


( 11 ) 


yields 
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( 12 ) 


m 


change ,j 


_ j^Wi 

dt 


dx 


u 


pj 


where u p } is the velocity of droplets in the j™ size class. This velocity is calculated 
through a drag analysis developed later in the report. 

Equation 12 can be converted to a change in mass flux by multiplying m change j by 
the total number of droplets in the j^ 1 size class passing through the control volume per 
unit time, N/Xj. Substituting this result, Eq. 9 becomes 


dny 

dx 





u 


P.J v 


dt 


/ I.H 


+ 


" dnyp 

v dt 


Al 


+ 




( 13 ) 


The value of N/t, for each droplet size class can be determined from the total initial 
slurry mass flux and normalized droplet size distribution. Given an initial normalized 
droplet size distribution, it is necessary to determine N/x total before N/Xj can be found for 
the individual size classes. N/x total can be found from the following mass conservation 
expression: 


m 


N 


sluny .initial 


Zf*i7P< 

total j=l \ ° 


(14) 


where Nj represents the percentage of total droplets in a given size class, given by the 
normalized size distribution. The variables, d J init and p, represent the initial diameter of 
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the jth size class and the initial slurry density, respectively. With N/x tota , known, N/Xj for 
each size class can be calculated from 


N 


x 


j 



total 


(15) 


Neglecting droplet breakup and secondary atomization, N/Xj is constant throughout the 
chamber This leaves the following three variables for each size class to be determined 
from the droplet gasification models: (dm*^ t /dt ) LH ’ (^drop./dt)^ . 

Energy Conservation : The steady-state energy balance for the system, including radiation 
heat losses, can be expressed as 


d(dih) gas = 

dx fr 

where h represents total specific enthalpy (h chem + h^) and q r is the radiation heat flux 
from the flow to the wall. Since gas flow optical properties are dominated by the small 
A1 2 0 3 fume particles 12 , the radiation term is independent of the individual droplet size 
classes and is not included in the summation term. 

Each of the bracketed terms on the right hand side of Eq. 16 can be expanded as 


d(mh) 


dx 


i-h , dfmh)^ 


dx 


d(mh) 

dx 


AKO, 


-q 


(16) 


d(mh) 

dx 


II 

. dh .dm 
m — + h — 

J 

dx dx _ 


( 17 ) 


Substituting for dm/dx yields 


12 


( 18 ) 


d(mh) 

dx 

J 


. dh f N 1 dm* "l 

m — + h L 

dx tu dt 


Because the enthalpy required to heat the liquid hydrocarbon from 300K to the 
hydrocarbon boiling temperature (438.9 kJ/kg at P = 6 atm) is comparable to the 
hydrocarbon latent heat of vaporization (~286.7 kJ/kg at P = 6 atm), hydrocarbon heatup 
must be accounted for in the hydrocarbon vaporization model developed later in this 
report. In the vaporization model, the bulk droplet temperature is assumed to remain at 
the initial droplet temperature. Heat transfer from the gas heats a thin surface layer of 
liquid hydrocarbon from the initial temperature to the boiling point and vaporizes this 
hydrocarbon layer. Using the constant bulk droplet temperature assumption yields a 
constant liquid hydrocarbon specific enthalpy, which reduces Eq. 1 8 for the liquid 
hydrocarbon to 


d(mh) LH 

— Vi 

N 1 

dx 

j 

t u p dt 


(19) 


where h LHTinit is the liquid hydrocarbon specific enthalpy at the initial droplet temperature. 

Following hydrocarbon burnout, an agglomerate of aluminum particles 
remains. 16 ’ 17 The aluminum agglomerate temperature rises from the hydrocarbon boiling 
temperature, through the aluminum melting point, to the aluminum boiling temperature as 
heat is transferred from the gas flow to the agglomerate The enthalpy required for this 
temperature increase is significant compared to the total system enthalpy and the 
aluminum enthalpy of vaporization, h^. Consequently, agglomerate ignition/heat up 
should be modeled. Presently, details of the agglomerate ignition/heat up process itself are 
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neglected but the enthalpy transfer required by the process is not. Upon hydrocarbon 
burnout, the aluminum agglomerate temperature is assumed to jump immediately to the 
aluminum boiling point, and the enthalpy required for this temperature jump is subtracted 
from the gas flow enthalpy. 

Approximating the agglomerate ignition/heat up process as a step change 
significantly underpredicts droplet combustion times and is intended only as a temporary 
measure until other aspects of the code have been developed. Following this 
development, an aluminum agglomerate heat-up model will be incorporated into the code 
After the step heat-up, Eq. 1 8 for the aluminum mass flux becomes 


d( mh)^ 
dx 


— h 

N I M 

“ "Al.Tbju 
J 

. x u p dt 


( 20 ) 


where h AI T is the liquid aluminum specific enthalpy at the boiling temperature. 

Since A1 2 0 3 agglomerate forms on the aluminum droplet surface, the A1 2 0 3 
agglomerate is maintained at the aluminum boiling temperature as long as any aluminum 
remains in the droplet. Ignoring the time period after aluminum burnout when A1 2 0 3 
agglomerate temperature can vary, Eq. 18 reduces to 


d(mh) 




dx 


- ^AIjOj.T^ 

N 1 

^^<fe-op,AJ 2 0 3 

T U 

dt 

J 

_ r 



( 21 ) 


for the A1 2 0 3 agglomerate, where h A]j0j is the A1 2 0 3 specific enthalpy at the aluminum 

boiling temperature. Substituting Eqns. 19, 20 and 21 into Eq. 16 yields the final system 
energy balance: 


14 


d(mh) g 

dx 


J=n N, 


dm 






drop.LH 

dt 


+ h 


dm 


AI.Tk, 


drop. A] 

"dt 


+ h 


dm, 




drop,AI 2 0 3 

dt 


-q r ■ (22) 


Radiation Heat Transfer : Radiation from the solid combustion products to the chamber 
wall is a participating medium phenomena requiring the solution of the radiative-transfer 
equation. Expressed in cylindrical coordinates, this equation is 


sinG 


cos<j> 


^I(r,9,<t>) sin <t» a(r,e,«t>) 


dr 


Aj) 


l(r,0,4)) = (l-aK,p(r) + ^ jll(r,0',4>')smO'd0'd<t.' , (23) 


where l(r,0,<j>) is the radiation intensity, r is radial distance from the centerline, and the 
variables, 0 and <J>, represent the angular spherical coordinates. I bp (r) is the blackbody 
radiation emitted by the fume particles, and Q 0 is the scattering albedo, defined as 


n 0 =-^— (24) 

a + o s 

where a and a s are the absorption and scattering coefficients, respectively. The absorption 
and scattering coefficients are currently unknown properties that must be determined 
before the transfer equation can be solved 

Given a scattering albedo, the transfer equation can be solved numerically, but 
requires a great deal of computer time. Fortunately, the transfer equation can be 
simplified and solved using one of several approximations based on the value of the optical 
thickness, K D , which is defined as 
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s s 

K d = J (a + a s )ds = J Kds , 

0 0 


( 25 ) 


where S is the characteristic radiation path length. IfK, defined as (a + cr 5 ), is constant 
along this path, K D can be expressed as 

K d = kS . (26) 

A1 2 0 3 smoke optical properties are used in evaluating the transfer equation since 
the smoke particles dominate the gas flow optical properties as discussed earlier in the 
report. Parry and Brewster 12 determined the optical properties of A1 2 0 3 smoke produced 
by a solid propellant containing 20% aluminum by mass, burning at a 1.8 MPa pressure. 
Optical thicknesses, K D = 1 .28 at X = 632.8 nm and K D = 1.21 at X - 1064 nm, were 
found for the 1 mm thick smoke region Assuming a constant k in the transmission 
direction, Eq. 26 produces k = 1280 and k = 1210 for Xj and Ao, respectively. 

Although the above k values were obtained from a solid propellant flame, they 
should provide a first approximation of the gas flow optical thickness in the model 
combustion chamber. Given that k is a function of fume volume fraction, k actually may 
be greater for aluminum slurry combustion than for the above solid propellant due to the 
higher slurry aluminum mass percentage. As a conservative estimate. Parry and 
Brewsters' k values will be used in estimating a minimum gas flow optical thickness. 

In the model, chamber radius, R = 0.025 m, is the characteristic path length Using 
this radius, the optical thickness is between 30 and 32. These values are much greater 
than K d = 2.0, which is considered the lower limit for an optically-thick medium. 18 
Therefore, as a first estimate, the gas flow can be treated as an optically-thick medium, in 
which local radiation heat transfer is only influenced by the immediate surroundings 
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Based on this optical thickness, a diffusion approximation 18 can be used to simplify the 
radiative transfer equation 

The cylindrical P, diffusion approximation 19 has been chosen since it provides 
greater accuracy than other approximations, yet is still simple to incorporate In the P, 
approximation, the radiant intensity is treated as a series expansion, 

l(r,e,+) = £ ix(0,*) <Pjr) , (27) 

n-0 

where T^r) are unknown functions, and Y^O,^) are defined as 


P„”(cose)e w (28) 

y(n + mj! 

In Eq. 28, P n m ( cos©) are Legendre polynomials and i = V-T. 

Eq. 27 is substituted into the transfer equation (Eq. 23) and the resulting equation 
is then mathematically manipulated, yielding the following expression: 




d^n+Um+l ( r ) , m+ l U J /,\ 

, + ^n+],m+l V r t 

dr r 


B„ 


d v P nt , m _, (r) m-1 


dr r 


+'* / ™ (O - [( 1 ■ - )I b (r ) + ^oo ]5o„5 


D„ 


d'W.fr) m -1 


dr 


dr 


'W.W 


*W,(r) 


Om 


= 0 , 


(29) 


where is the Kronecker delta. Since Y nm (0,<j>) does not equal zero, the expression in 
the { } brackets must equal zero for all n and m 
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Assuming that the radiant intensity is axisymmetric about the chamber centerline 
(I(r,0,<J)) = I(r,0,-<|>)), and that the radiant intensity does not vary along the chamber 
centerline (I(r,0,<j)) = I(r,7r-0,<j>)), Eq. 29 yields a set of N+l coupled equations where the 
number of equations is related to the order of the approximation; i.e. n=0,l,...,N In the 
Pi approximation, this series is terminated at N=l, resulting in the following two coupled 
equations: 


d ; f M (r) | 1 dTjr) 
dr 2 r dr 


3(l-fi 0 )'F„(r) = -3(l-n 0 )l l| ,(r) 


(30) 


and 


+ J 2T M (r) = 0 , (31) 

dr 

where I bp (r) is the blackbody intensity emitted by the fume particles. Since the combustor 
model is one dimensional, particle temperature is independent of r; therefore I b p (r) is 
independent of r 

The physical significance of 'E 00 (r) and T^^fr) can be found from the definitions of 
incident radiation, G(r), and net radiative heat flux, q r (r), respectively. The incident 
radiation is defined as 


2 n n 

G(r)= J Jl(r,0,(j))sin0d0d<t) 

(j)--oe^o 


( 32 ) 


Substituting the ?! series expansion for l(r,0,<j>) and integrating yields 
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( 33 ) 


«, _G(r) 

°’° ' 4k ' 

Similarly, the definition of net radiative heat flux, 


2 71 7T 

q r (r)= [ f l(r,0,<}>) sin 2 0cos<t»d0d<j> , (34) 

<J>=O0=O 


results in the following expression for 'P , 


(r) 



(35) 


Replacing T 0 0 and T, , in Eq 30 and 31 produces the following system of 
radiative transfer equations: 


^ll + ii£W_3(l_n„)G(r) = -12l(!-n„)l tp (36) 

dr' r dr 

^ill+3q r (r) = 0 . (37) 

dr 

Equation 36 can be solved for G(r) given the appropriate boundary conditions, and Eq 37 
can be solved for q r (r) using the solution for G(r). Using the fact that G(r) is axisymmetric 
about the centerline (r=0) and assuming that the chamber wall (r=R) behaves as a 
blackbody, the boundary conditions can be expressed as 
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(38a) 


dG(0) 
dr 

G ( R ) + l^) =4 , t I l .„ al (38b) 

3 dr 

Equation 36 can be transformed to a modified zero-order Bessel equation by using the 
variable substitution 


5 




where 




(39) 


Using this substitution and integrating Eq. 36 subject to the boundary conditions, the 
following expression for incident radiation is obtained: 


G(r) = 4 7cl bf + 





(40) 


I 0 (mr) and I,(mr) are the modified Bessel functions 


i„(mr) 


- (mr / 2 ) 21c+v 
^ 0 k!r(k + v+l) ’ 


(41) 


where v equals 0 or 1 and m = 1 /£. Substituting Eq. 40 into Eq. 37, the following 
expression for q r (r) is obtained: 
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( 42 ) 


4ro('C-' 

C)i, 

K) 

3$I 0 | 

[%) 

l + 2I,( 

>0] 


As mentioned above, the absorption and scattering coefficients, a a and cr s , are still 
unknown and must be determined prior to solving for q r (r). Consequently, radiation heat 
transfer is currently neglected in the system energy balance (Eq. 22) but will be included 
upon successful determination of a a and o s 


Momentum Conservation : The gas-phase momentum equation in the combustion chamber 
is trivial, assuming a negligible chamber pressure gradient and no body forces. The 
negligible chamber pressure gradient condition should be accurate for current conditions, 
but does require the assumption of no wall frictional losses or pressure drops due to flow 
acceleration. However, the momentum equations governing the hydrocarbon, aluminum, 
and agglomerate A1 2 0 3 mass fluxes are significant. In a given slurry droplet, there is no 
slip between the hydrocarbon and the aluminum before hydrocarbon burnout, and no slip 
between the aluminum and the A1 2 0 3 agglomerate afterwards; therefore the hydrocarbon, 
aluminum and A1 2 0 3 agglomerate momentum equations in a given droplet size class can be 
reduced to a single momentum equation. Virtual mass and Bassett forces can be neglected 
since the particle density is much greater than the gas density. The particle momentum 
equation for a droplet size class can be expressed as 

F = ma , (43) 
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where F is the drag force on a particle, m is the particle mass and a is the particle 
acceleration, du p /dt. Substituting for the drag force using a drag coefficient, C D , and using 
the chain rule to relate dup/t to dUp/dx results in the following form of Eq. 43: 

du p _ 3 PgasC D (u g ~ u P )|u g -u p [ (44) 

dx 4 P p d u p 

In this equation, p gas is the gas density, u g the gas velocity, u p the droplet velocity, p p the 
droplet density, and d the droplet diameter 

The drag coefficient is approximated as that of a sphere using the following 
correlation 20 : 


C 
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(45) 


where the Reynolds number. Re, is based on the slip velocity between the gas and the 
droplet. 

Slurry Combustion Overview : The physics of slurry droplet combustion and a probable 
secondary atomization process are described in detail in other research 6 > 71 61 7 21 22 In 
brief, when a slurry droplet is exposed to a hot ambient environment, heat transfer from 
the gas flow causes liquid hydrocarbon to vaporize and the droplet surface to regress. As 
liquid hydrocarbon vaporizes at the droplet surface, the aluminum particles that were 
suspended in the now vaporized hydrocarbon remain behind at the droplet surface, causing 
the number density of aluminum particles at the droplet surface to increase. 
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If the initial droplet is large enough, and there is a sufficient number of aluminum 
particles in the droplet, the number density of aluminum particles at the surface will 
increase until a rigid shell is formed by the particles coming into contact with each other. 
Some time interval after this rigid shell formation, surfactant pyrolysis causes the shell to 
become impermeable, and hydrocarbon vaporization to cease. Heat transfer to the slurry 
droplet continues, causing droplet internal pressure and rigid shell stresses to rise until the 
shell fails. Shell failure shatters the droplet, producing a number of smaller droplets which 
might repeat this secondary atomization process if the new droplets meet the minimum 
criteria for rigid shell formation 7 21 

After hydrocarbon burnout occurs in a droplet or droplet fragment, the remaining 
agglomerate of aluminum particles heats up, melts to form a single molten aluminum 
droplet, and bums. 

Secondary Atomization Modeling : Cho and Takahashi's shell formation model 21 is used to 
predict the droplet diameter at which rigid shell formation occurs for each droplet size 
class in the combustor model. When running the combustor code, it is assumed that 
secondary atomization occurs when the droplet diameter reaches the predicted rigid-shell 
diameter. Although secondary atomization actually occurs some time after rigid shell 
formation, the time interval presently is not known and is therefore neglected in the 
engine model. However, future experimental efforts should provide an estimate of this 
time interval, and then the time interval will be included in the engine model. 

Particle size distribution after secondary atomization is also currently an unknown 
to be determined from experimental measurements and is treated as a system variable. In 
the combustor model, a droplet undergoing secondary atomization is presently assumed to 
shatter into a specified number of equal-size secondary droplets. Defining the 
fragmentation ratio, P, as the number of secondary droplets produced per initial droplet, a 
new value of N/t: can be expressed as 
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( 46 ) 


N 

T 


j.new 



Knowing the slurry mass flux of a given size class, m sluny /r a new size class droplet 
diameter, dj new , can be found from the mass conservation expression. 


m 


N 


slurry j 


j.new 


-d 3 
6 J ' new 


( 47 ) 


Hydrocarbon Vaporization : Because of the close proximity of the slurry droplets to each 
other, the droplets are not surrounded by individual flames. Therefore, droplet 
combustion can be modeled as an evaporation process. Droplet heat up is approximated 
by assuming that heat transfer from the gas flow only affects a thin layer of liquid at the 
droplet surface. All of the heat transferred to the droplet from the gas flow is consumed in 
heating the thin surface layer from the initial droplet temperature to the liquid boiling point 
and vaporizing the layer. Consequently, the bulk droplet temperature remains at the initial 
droplet temperature. Assuming the following: 

• Quasi-steady droplet evaporation 

• Lewis number, Le = 1 

• Uniform droplet temperature equal to initial droplet temperature 

• Constant thermophysical properties, 

the gas-phase energy equation can be expressed as 
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(48) 


dr 


dX 


gas I _ 


dr 


: £t 

dr 


i y mC p,g^ 

where Z s 


4rck 


gas 


Integrating Eq. 48 using the boundary conditions, 

T (r — » x) = T a and T(rJ = T b , (49) 

where r s is the droplet surface radius and T h is the hydrocarbon boiling temperature, gives 
the gas-phase temperature distribution. 


x (T.-T b )exp(-Z/r)-T.exp(-Z/r,) + T b 
r l-exp(-Z/r s ) 


(50) 


The droplet surface energy balance is presented in Fig. 1 and can be expressed as 

m[h fg LH + J t *" c pXH dT) = q cond ( 5 

where q cond is the heat conducted from the gas-phase to the droplet, hf &LH is the 
hydrocarbon enthalpy of vaporization and c p LH is the liquid hydrocarbon specific heat 
Using Fourier's law for q cond , Eq. 51 becomes 



fg.LH 



= 4 ^^^ 



(52) 


Substituting the gas-phase temperature gradient, dT/dr, at the droplet surface and solving 
for m, yields 
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m = 


l^ln 


"p,gas 


c pi »( T .-Tj 


l>,».L, + [ I "V» dT 

* Lmi 


+ 1 


(53) 


Writing a mass balance for the droplet produces the following result: 


d m drop.LH _ ^Ttk^Ts ^ 

dt c„ _ 


c p ^(T.-T b ) 


^tg.LH + 


r 

JT™, 


+ 1 


(54) 


S^dT 


which is used in the system mass and energy governing equations. 

Aluminum Combustion : As mentioned previously, aluminum combustion proceeds 
through two different oxidation mechanisms; the first is aluminum vapor oxidation, and 
the second is droplet surface condensation/oxidation. This dual process phenomena 
requires a complex model for accurate results. Work is in progress to incorporate a 
detailed aluminum combustion model 17 - 23 in the combustor code. However, the aluminum 
combustion process is currently treated using a simple combustion model. In this model, 
all heat released by surface oxidation/condensation goes into the droplet with none lost to 
the gases surrounding the droplet. The following assumptions are also used: 

• Quasi-steady droplet evaporation 

• Uniform droplet temperature equal to aluminum boiling temperature 

• Lewis number, Le = 1 

• A1 2 0 3 surface agglomeration does not interfere with aluminum vaporization 

• Constant thermophysical properties 
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The gas-phase sensible energy balance is similar to that of the hydrocarbon droplet 
yielding a gas temperature distribution, 

(T f -Tjexp(-Z/r)-T,exp(-Z/r,) + T b 
l-exp(-Z/r s ) 

where Z and r s now apply to the aluminum droplet and T f is the aluminum flame 
temperature. 

However, the droplet surface energy balance is modified by the condensation of 
A1 2 0 3 on the surface. A schematic of the energy balance is presented in Fig. 2. 

Employing the assumption that all heat produced by the condensation/oxidation reaction 
goes into the droplet yields 

^Al^fg.AI = ^.gas + ^AljOjhfg.AUOj ’ (^ 6 ) 

where m^ 0 is the rate at which A1 2 0 3 condenses on the droplet surface. The A1 2 0 3 
condensation rate, 0 , can be expressed as a fraction, r|, of the aluminum vaporization 
rate, m AI . Substituting x\m M for m A]j0j and rearranging Eq. 56 produces 

®A1 (^fg.Al — ^l^tg.AljOj ) ^Ic.gas 

A value of rj = 0.4 was determined from previous aluminum combustion work involving 
both experimental and numerical analyses. 17 

Substituting Fourier's law for q c gas as was done for the hydrocarbon and solving 
for m A1 yields the following: 
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(58) 


47tk„ r 


m Ai =■ 


-In 


P.gas 


C p,gas( T f T b) | 

(hfg .aj ~ nh fe ,Ai,0 5 ) 


Mass conservation for the aluminum droplet can be expressed as 


dm ^p.AI _ 4 ^ k g, S r s ^ 

dt c_ 


C p,gas (Tr "^b ) 
(hfg.AI _ AIjOj) 


+ 1 


(59) 


It should be noted that h fg AUOj is not a true enthalpy of vaporization since A1 2 0 3 

does not exist in a gaseous state Rather, it dissociates into 2A1 + 3/20 2 . Therefore 
hfg Ai,o, is actually the enthalpy released by the chemical reaction, 2Al (g) + 3/20 2 ->Al 2 0 3(l) , 

occurring at the droplet surface temperature 


In both the hydrocarbon and the aluminum droplet models, the thermophysical 
properties, c p and k gas , are assumed to be constant However, both of these properties are 
actually functions of gas temperature, which varies with radius from the droplet surface to 
the ambient gas environment. Thus, it is necessary to determine an average value for c p 
and k gas . Following the work of Law and Williams 24 , the following values of c p and k gas 
are used: 


c„=<v(T) 

k„=0.4k I (T) + 0.6k.(T) 


(60) 


where T = (T b + T J/2 . In the above equations, the subscript f denotes the gaseous fuel 
property, either hydrocarbon or aluminum. 


30 


Solution Method 

Equations 13, 22, 44, 54, and 59 form the governing equation set and are 

numerically integrated in the axial direction using the IMSL Dverk integration routine 25 to 
determine m g , (mh) g , u p , m LH , m AI , and m AJ/)( These values are then used to 

determine the gas temperature, T g (x), the radiation heat flux, q(x), final AI 2 O 3 
agglomerate diameter, and chamber burnout length. Initial conditions are supplied at the 
injector face. The solution of these equations requires the knowledge of the following 
properties at each solution point: 

• Gas-phase T, p, p and k 

• Liquid hydrocarbon h, h^, p, and T b 

• Vapor-phase hydrocarbon k and c p 

• Liquid aluminum h, hfg, p, and T b 

• Vapor-phase aluminum k and c p 

• A1 2 0 3 agglomerate h at T b ^ and p 

Physical Properties 

Gas-phase temperature, density, and composition are calculated using the 
STANJAN chemical equilibrium subroutine . 26 Temperature-dependent gas-phase 
conductivity and viscosity are currently approximated as those of 0 2 using third-order 
polynomial curve fits of Svhela's tabulations . 13 Coefficients and valid temperature ranges 
are presented in Tables la and lb. 

The gas mixture is approximated as a single gas because of the complexities of 
calculating gas mixture properties. Gas species interactions must be accounted for to 
accurately determine these properties. In the future, gas-phase conductivity and viscosity 
will be calculated using a property subroutine that accounts for these mixture effects. 
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Table I a. 

0 7 Conductivity, k = A + BT + CT 2 + DT 3 (W/m-K) 


A 


9.6789 x 10" 4 
-9.8981 x 10* 3 
3.016 x 10' 2 
1.0799 x 10' 3 
6.791 x 10' 2 
-3.3408 x lO- 2 
1.4209 x 10‘ 1 


9.6973 x 10’ 5 -3.3806 x 10*° 1.072 x 10' 

1.2178 x 10 -4 -4.8 x 10‘ 8 1 0847 x 10 

3.888 x lO' 5 9.2148 x 10' 9 -2.3244 x 10 

8.4393 x 10‘ 5 -1.4793 x 10' 8 1.937 x 10* 

4.6007 x 10 -6 1.7085 x 10' 8 -2.3244 x 10 

1.0636 x lO" 4 -1.7201 x 10’ 8 1.54% x 10' 

4.0726 x 10' 5 2.4058 x 10' 8 -2.3244 x 10 


100 - 1000 
1000 - 1500 
1500 - 2000 
2000 - 2500 
2500 - 3000 
3000 - 3500 
3500 - 4000 



0? Viscosit} 

Table I b. 

, n = A + BT + CT 2 + DT 3 (kg/m-s)- 10 6 

A 

B 

c 

D 

Temp Range (K) 

4.5863 

6.0007 x 10‘ 2 

-2.0733 x 10' 5 

4.2703 x 10* 9 

300-2000 K 

16.601 

3.5227 x lO* 2 

-3.3703 x W 6 

2.2687 x 10* 10 

2000-4000 K 
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This property subroutine has already been developed for previous work in modeling single 
aluminum droplet combustion. 17 

The JP-10 properties (p, c p , h,- g , and T b ) are taken from Szekely. 27 Liquid JP-10 
density was determined as 


Pjp _ 10 = 1166.4 -0.792 T(k) kg/m 3 . (61) 


Similarly, liquid specific heat was correlated as 

c p =257.32 + 4. 5187-T(K) J/kg-K . (62) 


JP-10 boiling temperature was determined from the following curvefit of temperature as a 
function of vapor pressure: 


-4704.2 

ln(p(Pa) / 3.069 x 10 9 ) 


K 


(63) 


By setting the pressure equal to chamber pressure, Eq. 63 yields the boiling temperature. 
Presently, the enthalpy of vaporization as a function of chamber pressure is unknown. 
Therefore, the enthalpy of vaporization at 1 atm will be used until a more accurate 
enthalpy is determined. Using the Clausius-CIapeyron relationship, Szekely 27 found the 
JP-10 enthalpy of vaporization to be 


h fg,jr-io = 286 7 kJ / kg 


(64) 


The liquid aluminum properties (h, hjg, and T b ) are determined from the JANAF 
tables 28 , as developed below, and gaseous aluminum properties (c p and k) are curvefits of 
Svehla's tabulations 13 , presented in Tables Ila and lib. 
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The JANAF tables provide the specific enthalpy of both liquid and gaseous 
aluminum over a significant temperature range, permitting the aluminum enthalpy of 
vaporization, h fi , A] , to be calculated over the temperature range. The following third-order 
polynomial curve fit for hf &A] was then determined: 

h tg .Ai = 11 905x1 ° 6 " 400 591-2 4365x1 ° 3T " + 3 4091 x10 " 7 73 J/k S (65) 


and substituted into the Clausius-Clapeyron relationship. 


dP hf g 

dT Tv fg 


( 66 ) 


Assuming that the vapor specific volume, v g , is much larger than the liquid specific 
volume, Vj, and that the vapor behaves as an ideal gas, the Clausius-Clapeyron relationship 
can be expressed as 


= — (l 1.905 x 10 6 - 400. 59T- 2.4365 x 10' 3 T 2 +3.4091 x 1 0“ 7 T 3 ) , (67) 

RT- 

where P is saturation pressure, T is saturation temperature, and R is the specific gas 
constant. Rearranging and integrating yields the following expression: 


dP 

dT 


Rln 


vP \J 


= — 11.905x 10 ( 


T T 


-400.3SH In 




v T i y 


-2.4365x 10' 3 [T, -T,] + 1.7045 x 10‘ 7 [T 2 2 — T, 2 ] . (68) 
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Substituting the known saturation conditions, = 1 atm, and T sat = 2790.812 K, for P, 
and T, in Eq. 68 provides a relationship for the saturation pressure, P 2 , as a function of T 2 . 
However, Equation 68 cannot be rearranged to yield T 2 as a function of P 2 Therefore, T 2 
is determined through a Newton-Raphson solution of Eq. 68. If a given saturation 
pressure, P 2 , is equal to the combustion chamber pressure, then the corresponding 
temperature, T 2 , is the aluminum boiling temperature. Once T b A] has been found, h^ is 
determined from Eq . 65. 

The liquid aluminum specific enthalpy, h y, is determined through a third-order 
polynomial curvefit of JANAF data, resulting in the following expression: 

h A1 = -29467 + 1 1 77. 5T + 4. 5 1 06 x 1 O ' 5 T 2 - 5. 1 928 x 1 0" 9 T 3 (J/kg). (69) 

We assume the A1 2 0 3 agglomerate temperature is equal to the aluminum boiling 
temperature as long as aluminum is present in the droplet, after which the A1 2 0 3 
agglomerate is free to equilibrate with the gas flow temperature. Therefore, the A1 2 0 3 
agglomerate should exist only in the molten state in the combustion chamber since the 
A1 2 0 3 melting temperature (T m = 2315 K) is much less than both the aluminum boiling 
temperature (T b > 2790 K) and the gas temperature once the aluminum is consumed 
Following the work of Kirshenbaum and Cahill 29 , the molten A1 2 0 3 density, Pa]j03 , is 

modeled as: 


Pa ij 0 , = 5632 - 1. 127T(K) kg/m 3 . (70) 


The liquid A1 2 0 3 enthalpy is determined through a third-order polynomial curvefit 
of JANAF data: 


Vo 3 = 1 - 504d7 + 839 011 + 0 2335412 - 3 0298 xl0 " 5T ' ( J/k 8) ( 71 ) 
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Theoretical Results 


The one-dimensional engine code was exercised using the chamber diameter, 
pressure, flow rates, and aluminum mass loading presented in Table III These values were 
chosen to simulate Galecki's test conditions. 8 Presently, chamber length is varied to allow 
complete propellant combustion. 


Table m. 

Model Operating Conditions 

Chamber Diameter 

0.0522 m 


Chamber Pressure 

690 kPa 


Slurry Flow Rate 

0.00142 kg/s 


Aluminum Loading 

60% 


Oxidizer Flow Rate 

0.00312 kg/s 



A normalized droplet size distribution remains to be determined from other 
researchers' work on slurry atomization. A literature search has been started to find 
appropriate references on this topic. For the present time, the arbitrary normalized size 
distribution in Fig 3 is used 

Secondary atomization effects on combustion chamber gas and droplet velocities 
can be seen in Figs. 4 and 5. The data in Fig. 4 were calculated assuming no secondary 
atomization while the data in Fig. 5 assume a fragmentation ratio, P, of 40. In both cases, 
smaller droplets equilibrate more rapidly with the gas velocity than large droplets, and the 
gas velocity increases along the combustor axis due to decreasing gas density as 
temperature rises and mass addition to the gas flow from the slurry droplets. 

As expected, droplet velocities are the same in both cases until secondary 
atomization occurs, after which, the slope of the Fig. 5 velocity profiles decreases sharply 
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NORM. 



DROPLET DIAMETER (/xm) 


Figure 3. Normalized droplet size distribution used in combustor code. 
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VELOCITY (m/s) 



Figure 5. Gas and droplet velocities versus axial location for a 60 wt% A1 slurry 
assuming a secondary atomization fragmentation ratio of 40. 
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and the droplet velocities rapidly equilibrate with the gas velocity. This rapid 
equilibration is caused by secondary atomization shattering the initial droplets into small 
secondary droplets which equilibrate with the gas velocity more quickly than the initial 
droplet could 

Typical major and minor gas flow species mass fractions are presented in Figs. 6 
and 7, respectively. Again, a fragmentation ratio of 40 was used in calculating the data in 
both of these figures. The jagged shape of the mass fraction curves can be attributed to 
secondary atomization effects on burning rates and the sudden extraction of enthalpy from 
the gas flow to account for aluminum heat up. This second effect will be eliminated upon 
the addition of a proper aluminum heat up model 

From Fig. 6 , it is seen that the 0 2 mass fraction decreases steadily as droplet 
combustion progresses. However, the flow is oxidizer rich, resulting in excess 0 2 upon 
fuel burnout. C0 2 and H 2 0 mass fractions reach maximums at approximately 0.04-0.05 m 
and then decrease as increasing temperature causes dissociation, evidenced by increasing 
CO and OH mass fractions in Figs. 6 and 7. It should also be noted that the liquid AI 2 O 3 
mass fraction of 0.25 is much greater than the solid AI 2 O 3 mass fraction of 0.02. 

Gas temperature versus axial location for the same conditions is plotted in Fig. 8 . 
The maximum gas temperature of 3872 K is only slightly higher than the 3829 K predicted 
for equilibrium. This difference is attributed to the A 1203 oxide agglomerate temperature 
being fixed at the A1 boiling point instead of being permitted to equilibrate with the gas 
temperature Work is currently in progress to incorporate this equilibration process in the 
combustor model. 

In Fig. 9, chamber burnout distance is plotted versus fragmentation ratio, P, to 
illustrate the benefits of secondary atomization. It is readily apparent that only slight 
secondary atomization is required to yield significant decreases in droplet lifetimes. 

Higher secondary atomization intensities, represented by larger fragmentation ratios, have 
a lesser effect on burnout length because droplet lifetime is inversely proportional to 
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MASS FRACTION 



Figure 6. Mass fraction of major gas species versus axial location. Data is for a 60 wt% 
A1 slurry assuming a secondary atomization fragmentation ratio of 40. 
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MASS FRACTION 



X LOCATION (m) 


Figure 7. Mass fraction of minor gas species versus axial location for a 60 wt% A1 slurry 
assuming a secondary atomization fragmentation ratio of 40. 
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X LOCATION (m) 


Figure 8. Gas temperature versus axial location For a 60 wt% A1 slurry assuming a 
secondary atomization ratio of 40. 





CHAMBER BURNOUT LENGTH (m) 



FRAGMENTATION RATIO 


Figure 9. Chamber bumout distance as a function of secondary atomization 
fragmentation ratio for a 60 wt% A1 slurry. 
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droplet surface area, which increases as (J3) 2 ' 3 . It should be emphasized that the burnout 
distances presented in Fig. 9 neglect the A1 agglomerate heat up which is likely to increase 
these distances significantly . 

Final oxide agglomerate diameter as a function of secondary atomization intensity, 
P, is shown in Fig. 10. Similar to the trend seen in Fig. 9, small secondary atomization 
intensities significantly reduce final A1 2 0 3 agglomerate diameters, with greater atomization 
intensities providing decreasing marginal reductions in diameter However, final particle 
diameter is proportional to (1/P) 1/3 and not (1/P) 2 ' 3 as was droplet lifetime. 

The results presented in Figs. 4-10 demonstrate that the combustor model is 
performing as expected, given current model assumptions, and that only low secondary 
atomization intensities are required to significantly decrease chamber burnout length and 
final particle size. 


Experimental Efforts 

The main focus of our work during this reporting period was devoted to the 
theoretical model development described above. Experimental efforts consisted of 
improvements to the experimental apparatus described in earlier reports 30-33 . Previously, 
the transmitting and receiving optics were mounted on separate optical tables that were 
not isolated from the floor Building vibrations and people moving about the lab were 
sufficient to distort optical alignment. 

To correct this problem, both optical table tops and the burner were mounted on a 
single frame constructed of steel box beams. The entire assembly was then mounted on 
five isolation legs (Newport pneumatic isolation mount type XL-A). The optics were 
reassembled and the setup was recalibrated. Shakedown tests indicate that the system is 
correctly calibrated and that the apparatus is now sufficiently isolated from vibrations. 
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FINAL AGGLOMERATE DIAMETER (/im) 



Figure 10. Final oxide agglomerate diameter as a function of secondary atomization 
fragmentation ratio. Data is for a 60 wt% A1 slurry. 




Future Plans 


Plans for the next six months include the following: 

1 . Incorporate A1 heat up, gas mixture properties, radiation heat transfer, and improved 
aluminum combustion model in the one-dimensional rocket code. 

2. Determine whether existing rocket nozzle codes are suitable for modeling the 
expansion of slurry combustion products. If not, develop a one-dimensional nozzle 
code incorporating two-phase flow losses. 

3. Begin characterizing slurry and secondary atomization effects on engine 
performance 

4. Examine role of surfactants in secondary atomization by adding additional surfactants 
to existing slurries. 

5. Continue investigating the sample volume size probability distribution function and 
its sensitivity to particle index of refraction. 

6. Use existing experimental apparatus and additional imaging techniques to resolve 
post-secondary atomization droplet size distributions. 
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